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STABILITY OF FINITE DIFFERENCE SCHEMES 
FOR HYPERBOLIC INITIAL BOUNDARY VALUE PROBLEMS: 
NUMERICAL BOUNDARY LAYERS. 

BENJAMIN BOUTIN & JEAN-FRANgOIS COULOMBEL 


Abstract. In this article, we give a unified theory for constructing boundary layer expansions for dis¬ 
cretized transport equations with homogeneous Dirichlet boundary conditions. We exhibit a natural as¬ 
sumption on the discretization under which the numerical solution can be written approximately as a 
two-scale boundary layer expansion. In particular, this expansion yields discrete semigroup estimates that 
are compatible with the continuous semigroup estimates in the limit where the space and time steps tend 
to zero. The novelty of our approach is to cover numerical schemes with arbitrarily many time levels, while 
semigroup estimates were restricted, up to now, to numerical schemes with two time levels only. 


AMS classification: 65M12, 65M06, 65M20. 

Keywords: transport equations, numerical schemes, Dirichlet boundary condition, boundary layers, stability. 


Contents 


1. Introduction and main result 1 

1.1. Introduction 1 

1.2. Notations 2 

1.3. Main result 5 

2. Numerical boundary layers 6 

2.1. Formal derivation of the boundary layer expansion 6 

2.2. A preliminary result 7 

2.3. The leading boundary layer profile 10 

2.4. The first boundary layer corrector 11 

2.5. The approximate solution 12 

3. Proof of the main result 13 

3.1. The case of an incoming velocity 14 

3.2. The case of an outgoing velocity 15 

3.3. The semigroup estimate 18 

4. Example and counterexample 19 

4.1. A 4 time-step 5 point centered scheme 19 

4.2. The leap-frog scheme 21 

References 22 


1. Introduction and main result 

1.1. Introduction. The analysis of numerical boundary conditions for hyperbolic equations is a delicate 
subject for which several definitions of stability can be adopted. Any such definition relies on the choice of 
a given topology that is a discrete analogue of the norm of some functional space in which the underlying 
continuous problem is known to be well-posed. The stability theory for numerical boundary conditions 
developed in [GKS72], though rather natural in view of the results of [Kre70] for partial differential 
equations, may have suffered from its ’’technicality”. As Trefethen and Embree [TE05, chapter 34] 
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say: the term GKS-stable is quite complicated. This is a special definition of stability, that 

involves exponential decay factors with respect to time and other algebraic terms that remove it significantly 
from the more familiar notion of bounded norms of powers”. More precisely, the definition of stability in 
[GKS72] corresponds to norms of type for the numerical solution {t denotes time and x denotes the 
space variable), while in many problems of evolutionary type one is more used to the topology. In 

terms of operator theory, the definition of stability in [GKS72] corresponds to resolvent estimates, while 
the more familiar notion of bounded norms of powers corresponds to semigroup estimates. Hence a natural 
-though delicate- question in the theory of hyperbolic boundary value problems is to pass from GKS type 
(that is, resolvent) estimates to semigroup estimates. In the context of partial differential equations, this 
problem has received a somehow final answer in [Met 14], see references therein for historical comments 
on this problem. In the context of numerical schemes, the derivation of semigroup estimates is not as 
well understood as for partial differential equations. Semigroup estimates have been derived in [Wu95] 
for discrete scalar equations, and in [GGll] for systems of equations. However, the analysis in [Wu95] 
and [GGll] only deals with schemes with two time levels, and does not extend as such to schemes with 
three or more time levels (e.g., the leap-frog scheme). 

In this article, we focus on Dirichlet boundary conditions and derive semigroup estimates for a class 
of numerical schemes with arbitrarily many time levels. The reasons why we choose Dirichlet boundary 
conditions are twofold. First, these are the only boundary conditions for which, independently of the 
(stable) numerical scheme that is used for discretizing a scalar transport equation, stability in the sense 
of GKS is known to hold. The latter result dates back to [GT81] and is recalled later on. Second, 
homogeneous Dirichlet boundary conditions typically give rise to numerical boundary layers and therefore 
to an accurate description of the numerical solution by means of a two-scale expansion. We combine these 
two favorable aspects of the Dirichlet boundary conditions in our derivation of a semigroup estimate. 

The study of numerical boundary layers has received much attention in the past decades, including for 
nonlinear systems of conservation laws, see for instance [DL 88 , GS97, GHGOl]. As far as we know, all 
previous studies have considered numerical schemes with a three point stencil and two time levels. In 
this article, we focus on linear transport equations and exhibit a class of numerical schemes for which the 
homogeneous Dirichlet boundary conditions give rise to numerical boundary layers. The stencil can be 
arbitrarily wide. As follows from our criterion, the occurrence of boundary layers is not linked with the 
order of accuracy of the numerical scheme, which is a low frequency property, but rather with its high 
frequency behavior. For instance, the Lax-Wendroff discretization displays numerical boundary layers 
when combined with Dirichlet boundary conditions (and such layers have the same width as for the Lax- 
Friedrichs scheme) but the leap-frog scheme does not^, though both Lax-Wendroff and leap-frog schemes 
are formally of order 2 . 

1.2. Notations. We consider a one-dimensional scalar transport equation: 

(1.1) dtu + adxU = 0, t>0,x>0, 

where the velocity is a 7 ^ 0. The transport equation (1.1) is supplemented with an initial condition uq that 
belongs to a functional space that is made precise later on. In the case a > 0, that is, if we consider an 
incoming transport equation, we also supplement (1.1) with homogeneous Dirichlet boundary condition: 

( 1 . 2 ) u{0, t) = 0, t > 0 . 

The hnite difference scheme under consideration is assumed to be obtained by the so-called method of 
lines, see, e.g., [GK095]. In other words, we start with (1.1) and first use a space discretization. The 
latter is supposed to be linear with r points on the left and p points on the right. In other words, we 
consider some coefficients a-r,... ,ap, where p,r are fixed nonnegative integers, together with a space 
step Ax > 0, and approximate (1.1) by the system of ordinary differential equations: 

1 ^ 

'^3 + ^ ~ 0 ’ 

l=—r 

^The leap-frog scheme rather generates incoming highly oscillating wave packets, as explained at the end of this article. 


(1.3) 
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where Uj{t) represents an approximation of the solution u to ( 1 . 1 ) in the neighborhood of the point 
Xj := j Ax. The integers r,p are fixed by assuming a_r 7 ^ 0 and Up 7 ^ 0. The latter system of ordi¬ 
nary differential equations is then approximated by means of a (possibly multistep) explicit numerical 
integration method. We refer to [HNW93, HW96] for an extensive study of numerical methods for or¬ 
dinary differential equations. Applying a linear explicit multistep method to (1.3) yields the numerical 
approximation 

k k—1 p 

(1.4) a. + = 0 . 

(T=0 cr=0 i=—r 

with k > 1 and fixed constants ao,..., ak, /3o,..., f3k-i- The multistep integration method is normalized 
by assuming |ao| + |/5o| > 0 and au = 1. In (1.4), we have made use of the notation A := At/Ax for 
the so-called Courant-Friedrichs-Lewy parameter. In what follows, the parameter A is kept fixed^, and 
we consider the space and time grid xj := j Ax, t” := nAt for j,n G N. For notational convenience, we 
introduce the (dimensionless) constant r > 0 that satisfies 

(1.5) Ax = T |a| At. 

We keep At G (0,1] as the only small parameter and Ax G (0,1/A] varies accordingly. 

Since we are approximating the transport equation (1.1) on the half-line M"*", the space grid is indexed 
by N. This means that the numerical approximation (1.4) takes place for j > r. We then supplement 

(1.4) with homogeneous Dirichlet boundary conditions on the ’’numerical” boundary: 

( 1 . 6 ) = 0 , 0 < j < r — 1 , n>k, 

independently of the sign of a. The scheme (1.4), (1.6) is ignited by k initial data, which correspond to 
the approximation of the solution to (1.1) at times t^,, t^~^. For simplicity, we assume that the initial 
data for (1.4), (1.6) are given by the standard piecewise constant approximation of the exact solution to 
(1.1). In other words, we set: 

1 r^j+i 

(1.7) Uj / uo(x —at”)dx, j>0, n = 0,...,k — l, 

J Xj 

where the initial condition uq for ( 1 . 1 ) has been extended by zero to in the case a > 0 . 

The following two assumptions are the minimal consistency requirements for the numerical scheme 

(1.4) . 

Assumption 1.1 (Consistency of the space discretization). The coefficients a-r,... ,ap in (1.4) satisfy 


p 


(1.8) 

, 


i=—r 


P 

(1.9) 

£ai = a 

i=—r 


Assumption 1.2 (Consistency of the linear multistep integration method). The coefficients oq, ..., ak, 
Po,..., Pk-i of the time integration method in (1.4) satisfy 

k k k—1 

'^aa = 0, = 

a=0 a=0 cr=0 

In the case k = 1, that is for numerical schemes with two time levels, the normalization gives oq = 
Oil = Po = 1, and (1.4) reduces to the standard form 

p 

-Uj + X ai = 0 . 

l=—r 

^This assumption could be weakened by assuming that the ratio |a| At/Ax is bounded from below and from above, but 
we shall restrict to the more common case where the ratio is fixed for simplicity. 
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If p = r = 1, we obtain the class of three point schemes that encompasses both the Lax-Friedrichs and 
Lax-Wendroff scheme. 

As a direct consequence of the first consistency condition (1.8), it appears that the scheme (1.4) admits 
a conservative form in the following sense. There exists a linear numerical flux function F with real 
coefficients: 

p-i 

F{vj, . . . , ■= ^ ^ f1 

i=—r 


such that 


V 

(1.10) ^ ^ '^j+l — P{Fj—r+l: • • ■ ) '^j+p) ■ • • ) l) • 

i=—r 

In particular, (1.4) also takes the conservative form 
k k—1 

(1.11) + •' E'’- ('=’(““"+1. ■ ■ •. “SO - ■ ■ ■. “SOi)) = 0 ■ 

(T=0 cr=0 


From the second consistency condition (1.9), it follows that F{u,... ,u) = au for any n G M. This is the 
usual consistency property of F with the exact flux (u i—>■ au) of the transport equation (1.1) written as 
a conservation law. 


Our final assumption is the standard £^-stability assumption for (1.4) when the scheme is considered 
on the whole real line J G Z: 


Assumption 1.3 (Stability for the Cauchy problem). There exists a constant C > 0 such that, for all 
At G (0,1], the solution to 

k k—l p 

'y ^ Oia ^ Pa 'y ^ UGN, 

(T=0 f7=0 i=—r 


satisties 


k-l 

sup ^ Ax < C ^ ^ Ax \uj\‘^ 
j'ez cr=o jez 


As is well-known, assumption 1.3 can be rephrased thanks to Fourier analysis, 
introduce the function A defined by: 


More precisely, if we 


( 1 . 12 ) 


p 

Vz G C \ {0}, A(z) = ^ ai , 

i= — T 


then applying the Fourier transform to (1.4) yields for all ^ G M: 

k k—l 

^ a, jd^(0 + A PaA { F ^^^) j?^(0 = 0 , 

(T = 0 (7 = 0 

where u"' is the piecewise constant function that takes the value u” on the cell [j Ax, {j + 1) Ax). The 
stability assumption 1.3 is equivalent to requiring that there exists a constant C > 0 such that for all 
r/ G M, and for all given xq, ... , x^-i G C, the solution (xCT)CTeN to the recurrence relation 

k k—l 

VuGN, ag Xn+a + A /3a ^(e^’^) Xn+a = 0 , 

(7=0 (7=0 


sup |XnP < C ( |xo|^ H-h |xfc_i 

nSN ^ 


satisfies 
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In particular, the closed curve {—Ayi(e*^), r/ G M} should be contained in the so-called stability region 
of the numerical integration method, see [HW96, Definition V.1.1]. Observe now that the consistency 
assumption 1.1 introduced above can be rewritten under the form: 

(1.13) yi(l) = 0 and yi'(l) = a7^0. 

Since A vanishes at 1, 0 should belong to the stability region of the numerical integration method, which 
implies (see [HNW93, Chapter III.3]): 

k 

(1.14) y^aa^ ^0. 

a=0 

Remark 1.4. In the ease k = 1, the stability assumption 1.3 is equivalent to: 

(1.15) VzGSS |l-Ayi(z)| < 1. 

In particular, assumption 1.3 constraints the CFL number A to be ’’small enough”, and A{z) can not be 
a negative number. 

1.3. Main result. The main result of this paper is the following theorem. 

Theorem 1.5 (Semigroup estimate). Consider a linear seheme of the form (1.4) satisfying the consistency 
assumptions 1.1 and 1.2, the stability assumption 1.3 and the ’’dissipative” assumption 2.1 introduced later 
on. Consider an initial condition uq G //^(M+) for (1.1) such that 

fuo(O) = 0 , if a < 0 , 

l^uo(O) = Mo(0) = 0 , if a > 0 . 

Let T > 0 and, for At G (0,1], let us define Nt as the largest integer such that AtNx < T. Let also 
/i G [0,1/3]. Then there exists a constant C > 0, that is independent of T, At, fijUo such that the solution 
('“j)i>o,n>o to (1.4)-(1.6)-(1.7) satisfies 

(1.16) sup < c(||Mo|li2 (R+) + ||uo|Ih 2 (r+)) . 

n<NT jyQ 

Let us observe that (1.16) is compatible with the ’’continuous” estimate 

sup ||tt(t)||^2(R+) < C* ||rto||2,2(R+) , 

as At tends to zero. The role of assumption 2.1 is to derive a boundary layer expansion for (i^j )j>o,n>0) 
that is to decompose («”) as in [DL88, GS97, CHGOl] under the form 

where the boundary layer profile depends on the ” fast” variable j = Xj /Ax and has exponential decay 
at infinity, while the interior profile depends on the ” slow” variable xj. As follows from the analysis 
below, the derivation of such two-scale expansions is not linked to any viscous behavior of (1.4) (as the 
scaling Xj /Ax might suggest at first glance). 

The parameter p can diminish the T-dependence of the constants in (1.16). In particular, given any 
e > 0 and T > 0, there holds 

sup ^Axju^p <C'(||uo||i2(R+) + 2Ati-=||uo||^2(R+)) , 

n<NT ^.>0 ^ ^ 

for At sufficiently small (depending on T). 

Section 2 is devoted to the construction of boundary layer expansions for solutions to (1.4), (1.6). 
Theorem 1.5 is proved in Section 3 by means of a careful error analysis. We discuss some examples in 
Section 4 together with the relevance of assumption 2.1. 
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2. Numerical boundary layers 


2.1. Formal derivation of the boundary layer expansion. Our first goal is to understand when the 
numerical solution {u'^)j>o,n>o of the scheme (1.4), (1.6) can be approximated by an asymptotic boundary 
layer expansion: 

+ u^\j, n, j>0,n>0. 

In the latter decomposition, we expect to have fast decay at infinity. The functions and u'’* are 
to be dehned in such a way that represents an accurate approximation of (uj) as At tends to 0. 

Roughly speaking, the term takes care of the interior behavior of the solution far from the boundary, 
and involves the boundary layer correction that is localized in a neighborhood of x = 0 and matches 
the boundary conditions (1.6). 

We shall force the approximate solution to satisfy the initial conditions: 

(2.1) uf^ = u], j>0,n = 0, ...,/c-l. 

In this way, the error (uj’(f — u”) will satisfy a recurrence relation of the form (1.4), (1.6) with ’’small” 
source terms but will have zero initial data. We also expect the approximate solution to satisfy (1.6), or 
rather 


(2.2) ~0, 0<j<r-l, n>k, 

where, by ~ 0, we mean for instance that should be 0(At) on the boundary. 

For technical reasons that will be made precise in Section 3, we shall define the boundary layer term 
through a two term expansion of the form: 

involving a zero order term plus a first order corrector that will be used to remove part of the 
consistency error. 


We follow the discussions in [DL88, GS97, CHGOl] and briefly present hereafter a schematic derivation 
of the equations that will govern the three sequences n™*, and To that aim, let us introduce 

the following consistency error: 




At 


E 

V(T=0 


k-l 


0(7 U 


app 

'j,n+a 


+ A ^ 




a=0 


with j > r, and n > 0. 

• At a fixed positive distance from the boundary, the limit At —>■ 0 corresponds to j —>■ +oo and 
the boundary layer term becomes negligible with respect to The above consistency error 
reads (up to smaller terms) 

^ / fc fc-i p 

ej,n+k ^ ^ Z] + 

\(T=0 (T=0 i=—r 



This quantity will be of order 0(At) provided that u™* is a smooth solution to the continuous 
equation (1.1) (recall the consistency assumptions of the numerical scheme (1-4)). 


• Glose to the boundary, that is for a fixed index j > r, the limit At —?■ 0 makes xj tend to 
zero. If the interior solution is smooth enough, we get (recall that j is fixed) t"") = 

M“*(0,t"’) + 0{At) and u™*(0, = u“*(0,t"’) + 0{At). Then the consistency error reads^ (up 

to 0(1) terms): 

/ k k-l p 

£j,n+k - ^ ( X] + A ^ X] + £, r+'^) 

\cr=0 a=0 i=—r 


^Here we use the consistency conditions for the coefficients in (1.4). 
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Due to the consistency of the numerical integration method, and assuming that depends 
smoothly enough on the time variable, we get 


^j,n+k — 


1 

Ax 


'k-l 




(2.3) 


Vcr=0 / £=—r 

The first boundary layer profile needs therefore to satisfy the recurrence relation^: 

p 

^ + £,r) = 0, j>r,n>k. 

i=—r 

In terms of flux quantities, the relation (2.3) corresponds to requiring 


p-i 

E 

£=—r 


/,rr“’0(j+£ + r,r) = C 


n\ _ /^st 


(2.4) 


with an integration constant that only depends on n, but not on j. The constant is easily seen 
to be zero due to the required behavior of the boundary layer profiles at infinity. In addition, 
the boundary condition (2.2) imposes (up to an 0{At) term) the trace of on the numerical 
boundary: 

^n) ^ ^ 0 < j < r - 1 , n > 0 . 

We still keep the index j fixed and expand the consistency error at the following order with respect 
to At. Assuming that is smooth enough so that its associated consistency error is 0(At) up 
to the boundary, the overall consistency error reads (up to 0(At) terms): 


£j,n+k - i X] 1 


^k-1 


At 


a=0 


\cr=0 


i= — r 


We then require the first boundary layer corrector satisfy: 


(2.5) 


i=—r 


/k—1 

(i+<.n + T 


-1 


^a.u“-0(j,t-+") = 0, j>r. 


( 2 . 6 ) 


Vcr=0 / (T=0 

Since our analysis considers numerical schemes of order 1 or higher, the precise value of on the 
numerical boundary does little matter since any other choice than the one below will introduce a 
new 0{At) error that will just have the same order as the interior consistency error. For simplicity, 
we therefore require satisfy; 

u“d(j^ e) = 0^ 0<j<r-l,n>0. 

The above formal derivation of the profile equations (2.3) and (2.5) motivates the analysis of the 
recurrence relation (2.3). More precisely, we are going to determine the solutions to (2.3) that tend to 
zero at infinity. The precise definition of the approximate solution is given in subsection 2.5. 

2.2. A preliminary result. Let us recall that the function A, which is linked to the amplification matrix 
for the scheme (1.4), is defined in (1.12). The consistency assumption 1.1 implies that 1 is a simple root 
of A. The following assumption will turn out to be crucial in the forthcoming analysis. 

Assumption 2.1. The value z = 1 is the unique root of A on S^.' 

V0 G [—TT, tt] \ {0} , A(e*®)7^0. 

Remark 2.2. In the case k = 1, assumption 2.1 is obviously satisfied for every dissipative scheme (for 
which we recall that there exist c > 0 and A: G N* such that for all \d\ < vr, |1 — AA(e®®)| < 1 — 
However, we underline at this level that some non-dissipative schemes satisfy assumption 2.1 too, e.g. the 
Lax-Friedrichs scheme (that is considered in [CHGOlJj for which A(e*®) = cos0 — 1 — z Ao sin0j. 


^Recall that by our consistency and stability assumptions, the sum of the fia- is nonzero. 
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The main result of this subsection is the following Lemma. 


Lemma 2.3. Under Assumptions 1.1, 1.2, 1.3 and 2.1, the equation A{z) = 1 admits exactly R roots 
(with multiplicity) in D \ {0} = {z G C, 0 < |2| < 1} where 


R = 



if a < 0 , 
if a > 0 . 


Let us observe that in the case a > 0, r can not be zero and therefore one gets a nonnegative integer for 
R. Indeed, the value r = 0 is prohibited by the fact that the numerical dependence domain would not 
include the ’’continuous” dependence domain, see [CFL28]. 

The proof of Lemma 2.3 makes use of the following simple observation which we have not found in 
[HW96] and therefore state here. We keep the notations of [HNW93, Chapter III.2]. 


Lemma 2.4. Consider an ordinary differential equation of the form ij = f{y), and the explicit linear 
multistep integration method: 

k k—1 

(2.7) E Q((T Pn+cr — AR, 'y ^ fdfj fn+a ; 

a=0 tj=0 

with the normalization ak = 1, |q(o| + |/3o| > 0. Assume that the method is stable (in the sense of [HNW93, 
Definition III.3.2] j and that it is of order 1 or higher. Then the stability region for this method contains 
no positive real number. 


Proof. Following [HNW93, HW96], we introduce the polynomials 

k k—1 

q{X):=Y,»jX^ , ct(X):= 

j=0 j=0 

The assumptions of Lemma 2.4 can be rephrased as; 

^( 1 ) = 0 , g\l) = ail)^0, 

and g has no root of z satisfying \z\ > 1. In particular, must be positive for otherwise (recall ak = 1) 
g would have a real root in the open interval (1, +oo). We therefore have (t(1) > 0. 

For any given fJ. > 0, the real polynomial: 

P^{X) := g{X)-fia{X), 

has degree k and is unitary. It tends to +oo at +oo and T(i(l) = < 0. Hence vanishes in the 

open interval (1, Too) and p, does not belong to the stability region of the numerical method. □ 


Lemma 2.4 is consistent with the plots in [HW96] of the stability regions for the explicit Adams and 
Nystrom methods. Observe however that some stability regions may contain complex numbers of positive 
real part, e. g., the explicit Adams method of order 3. 

Proof of Lemma 2.3. Under assumption 2.1, A has no other zero on than z = 1 (with multiplicity 1). 
On the other hand, A admits a unique pole over C, at z = 0 and of order r (because we have a-r / 0). 
The cornerstone of the forthcoming proof is the residue theorem for meromorphic functions. Being given 
F a direct closed complex contour encircling the origin once and on which A does not vanish, then 

1 f A'(z) 

(2.8) - / ^ . ! dz = #{zeros inside F} — #{poles inside F} , 

ziTT Jy A[z) 

where zeros and poles are counted with multiplicity. The second integer on the right hand side equals r, 
and we intend now to compute Ry '.= ^{zeros inside F} thanks to an appropriate choice for the contour F 
(for which Ry = R). 
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The contour r^. Let us consider some parameter e G (0,7r/4] sufficiently small (to be determined later 
on), and let us define the contour F;. as but for a small chord avoiding 1, see Figure 2.1. More precisely, 
we consider the path F^ as the union F^p U with: 

Fe,! := j e*® , 0 G [e, 2 71 — e] > , '■= \ cose + iw , w G [— sine, sine] > . 



Figure 2.1. The integration contour Fg. 


Choice of the parameter e. Let us first observe that 1 is a simple zero of A so we can choose eo > 0 small 
enough such that, for any e G (0, eo], the number of zeros of A inside Fg equals the number of zeros of A 
in D \ {0}. 

Our goal now is to show that, for any sufficiently small e > 0, there holds 

> 0, if a<0, 

|±9yi(e±*^) > 0, if a>0, 

and for all z G Fg^ 2 ) aA{z) 0 M^. These properties follow from studying the variation of the function 
^^[(cose + ico). Namely, we compute 

— 9'yi(cose + itd) = 51?yi'(cose + fw) = a + 5R(yi'(cose + ico) — ^^(1)) • 
dij 

Consequently, if we assume a > 0, then (uj i—)■ ^^[(cose + iuj)) is increasing on [—sine,sine], while if 
we assume a < 0, then (a; i—>■ 9'yi(cose + i w)) is decreasing on [— sine, sine]. In any of these two cases, 
yi(cose + iio) is real for at most one value of to. 

We now observe that gd(cose) is real. In particular, for any 0 < jw] < sine, yi(cose + ico) belongs 
to C \ M and the sign property for is proved. Furthermore, using A'{1) = a, we find that 

>l(cose) is positive if a is negative while yi(cose) is negative if a is positive. We thus have, provided that 
e is sufficiently small, aA{z) 0 M"*“ for any z G Fg^ 2 - From now on, e is fixed and the latter properties hold. 


Application of the residue theorem. We denote hereafter log_ the principal complex logarithm with the 
usual branch cut along M_, and log_|_ the complex logarithm with a branch cut along M_|_. 

For any 2 ; G Fgp, one has A(z) / 0 (thanks to assumption 2.1) and A(z) 0 (for otherwise the 
stability region of (2.7) would contain —XA{z) G M"*"’*, which can not hold by Lemma 2.4). We can thus 
use log_ for computing the integral along Fg^i, and we get 


^ 

2iw Vr.,, Jf-a) 


;^(log_dl(e-‘') 

ZlTT \ 


log_yi(e*'')) . 


The integral along Fg ^2 depends on the sign of a. 
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Suppose o < 0. Then we know that for all 2 : G r£^ 2 ) •A{z) does not belong to 
use the log_ logarithm and derive 


1 


2 i TT 


- £,2 


A{z) 


dz = 


1 


2 i TT 


log_ yi(e*^) — log_ yi(e 


dz = 0, 


1 


2 i TT 


- £,2 


A(z) 


Summing the two contributions over Tg^i and re^ 2 ) we obtain 


R — r = 


1 


2iTT 


log+yi(e*^) - log_yi(e*^)) - 


1 


2 i TT 


We can again 


Summing the two contributions over T^^i and re ^2 we finally obtain 

1 f 

2iTr A(z) 

and R = r. 

Suppose now a > 0. Then we know that for all 2 ; G re^ 2 , does not belong to 
log_|_ logarithm and derive 

^ - log+ ■ 


We use the 


log_,_yi(e —log_yi(e 


The difference log_,_ — log_ equals 0 on 12+ := {2 G C, 9 2 > O} and equals 2i7r on 12_ := {2 G 
C, 92 < 0 }. To complete the proof, we recall that belong to f2±, and we thus get 

R — r = —1. 

□ 


2.3. The leading boundary layer profile. Using the flux function F, we can rewrite the boundary 
layer profile equations (2.3), (2.4), and introduce the following definition. 

Definition 2.5. Being given a real number u, we call a boundary layer profile associated with u 

a sequence that satisfies the following requirements: 

(i) Vo = ■■■ = Vr-l = -U, 

(a) F{u + Vj,... ,u + Vj+p+r-i) = F{u,... ,u), for all j > 0, 

(Hi) limj_,.oo Uj = 0. 


Let us comment some facts. The first point (i) above is related to the Dirichlet condition (2.4) with u in 
place of tt™*(0, f") (here the time variable is frozen). As a consequence of the linearity of the numerical 
flux F, the above condition (ii) reads 


which is equivalent to 
(2.9) 


p-i 

V j > 0, fe Vj+i+r = 0 , 

i=—r 

P 

V J > 0 , ae Vj+e+r = 0, 

i=—r 


if condition (Hi) is satisfied. Boundary layer profiles are therefore the zero-limit solutions to the linear 
recurrence relation (2.9) for which the r first terms of the sequence coincide. 


Definition 2.6. The set of all the values u such that a stable boundary layer associated to u exists is 
denoted 

Cnum = |u G M, G boundary layer profile associated with u| . 

This definition is the same as in [DL88, GS97, CHGOl]. The set Cnum encodes the so-called residual 
boundary conditions for (1.1) coming from the continuous limit At —)• 0 in (1.4), (1.6). We are now ready 
to prove the following result that characterizes the boundary layer profiles for the scheme (1.4). 

Proposition 2.7. Under Assumptions 1.1, 1.2, 1.3 and 2.1, there holds: 












NUMERICAL BOUNDARY LAYERS 


11 


• if a > 0, then Cnum = {0} and the unique boundary layer profile associated with 0 is the zero 
sequence (vj = 0 for all j > 0); 

• if a < 0, then Cn um = M and for any it G M there is a unique boundary layer profile 
associated with u, that decreases exponentially fast at infinity. We may write 

(2.10) Vj = uwj , j > 0 , 

where {wj)j^fq denotes the boundary layer profile associated with u = l. 

Proof. As explained above, our goal is to determine the zero-limit solutions to the recurrence (2.9) that 
satisfy condition (i) m. Definition 2.5. We thus look for the (stable) roots to the polynomial equation 

^ = 0 . 

l=—r 


Since this polynomial does not vanish at zero, its roots in D coincide (with equal multiplicity) with the 
zeros of A in D \ {0}. Lemma 2.3 precisely gives the number of such zeros. 

The zero-limit solutions to the linear recurrence (2.9) are spanned by the sequences (m = 1,..., r 
if a < 0 and m = l,...,r — lifa>0): 


( 2 . 11 ) 


{f zl)j(zn , 0<u<^j,l<i<g, 


where zi,... ,Zq denote the pairwise distinct zeros of A in D \ {0} and pLi,..., pLq their corresponding 
multiplicity. 

• We first assume a > 0. The subspace of zero-limit solutions to the linear recurrence (2.9) has 
dimension r — 1. Let n G M. We are looking for a sequence v = oJm such that 

Vo = ■■■ = Vr-i = —u, which is equivalent to 


/4‘’ 


^h-i) \ 


V 7(1) 7(^-1) 1 / 

—1 • • • ^r—l U J 


= 0 . 


The involved matrix in Mr,r(C) is invertible and therefore u = 0, u = 0. 

• We now assume a < 0. The subspace of zero-limit solutions to the linear recurrence (2.9) has 
dimension r. Let u G M. We are looking for a sequence v = Ylm=i such that vq = ■ ■ ■ = 

Vr-i = —u, which is equivalent to 


(Z, 


( 1 ) 

0 


7{r) • 




= —u 


Vzii\ ... zi"_V 






The involved matrix of Mr,r(C) is invertible and thus, for each given u G M there is a unique 
solution {uji,... ,L0r) G which determines the boundary layer profile associated with u. By 
linearity, this profile takes the form (2.10) and it is exponentially decreasing. 

□ 


2.4. The first boundary layer corrector. Our goal in this subsection is to construct a solution to the 
first boundary layer corrector equations (2.5), (2.6). In what follows, the function u*’*’*^ will be a boundary 
layer profile associated with some discretized trace of the exact solution to (1.1). In the case a > 0, there 
is no boundary layer profile but zero and the solution to (2.5), (2.6) is also zero. In the case a < 0, 
the space of boundary layer profiles is spanned by the sequence (rCj)jgN, and it is therefore sufficient to 
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construct a sequence that satisfies 

p 

( 2 . 12 ) yj>r, aiWj+i + Wj = 0, 

l=—r 

wq = ■ ■ ■ = Wr -1 = 0 , lim Wj = 0 . 

j-^oo 


Lemma 2.8. Under the assumptions of Proposition 2.1, in the case a <0, there exists a unique solution 
('uij)jgN to (2.12) and this solution decays exponentially fast at infinity. 

Proof. Uniqueness easily follows from the linearity of (2.12) and Proposition 2.7. As far as existence is 
concerned, we keep the notation of Proposition 2.7 and decompose the sequence {wj)j£n as 

r 

m=l 

where the are given by (2.11). Due to the linearity of (2.12), we first construct a zero-limit solution 

to the recurrence 

Vj>r, 

i=—r 

which is done by choosing of the form 

«''"■* = E fw"*" 4. 

/i=0 

and by identifying the coefficients go,..., (this procedure gives an invertible upper triangular system). 
Summing finitely many such sequences we get a sequence W that decays exponentially at infinity 

and that is a solution to the recurrence relation 

p 

V j > r , ^ + wj = t). 

l=—r 


The sequence (lUj) is obtained by correcting the initial conditions for (Wj), that is by choosing 


with 


r 

w:=W+Y,^m , 

m=l 



□ 


2.5. The approximate solution. Let us recall that the solution to (1.1), supplemented with the homo¬ 
geneous Dirichlet condition (1.2) in the case a > 0, is given by 

(2.13) u^^{x,t) = uo{x — at), x > 0, t>0, 

where the initial condition uq has been extended by 0 to in the case a > 0. This suggests defining the 
interior numerical solution as 

uo{x - at^) dx , j > 0 , n > 0 . 

In particular, (1.7) gives 

Vn = 0,.. ., A: - 1, V j > 0 , = u" . 
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In the case a > 0, there is no boundary layer and we define the approximate solution u^pp to (1.4), 
(1.6) as 

In the case a < 0, there exists a one-dimensional space of boundary layer profiles and we can also 
construct boundary layer correctors. In view of (2.4), we first need to approximate the trace of the exact 
solution and therefore set 

1 

Vn>0, '^n ■=-^ J uo{-at)dt. 

We now define the leading order boundary layer profile and first order boundary layer corrector 
as follows: 

bi,o jo, i>0, n = 0, ...,A:-1, 

U ■ \ 

' [UnWj, j > 0, n>k, 


bl,l 

^j,n ■ = 


0 , 


j >0, n = 0, ...,/c-l, 


ELo (ELo Wj, i > 0 , n>k. 

The approximate solution u'^pp to (1.4), (1.6) is then defined by: 


(2.14) 


app mt I Di,U I A Dl,l • \ n \ n 

^j n ■= Uj,n + Uj n + Ax , J > 0 , U > 0 . 


Thanks to our choice for the initial data, we again have: 


J A 0 , n = 0,..., A: - 1. 


3. Proof of the main result 

The error analysis uses the expression of the approximate solution {u^^n)j>o,n>o introduced in subsec¬ 
tion 2.5. We thus focus on the interior consistency error that is defined by: 


^j,n+k .— 


Al S 

Vcr=0 


fc-1 


app 

Uj n+a 


+ ^ /^o- ^ Rr 

(T=0 t=—r 


U 


app 

'j-\-^^n+a j ’ 


with j > r and n > 0, and on the boundary errors: 


Vj,n ■■= 0 < j <r - I, n>k. 

We recall that the approximate solution has the same initial data as the exact numerical solution 
(whatever the sign of a): 

EE = E ’ i A 0 , n = 0,..., /c - 1. 

Consequently, the error: 

ej-n:=EE-E’ JAO, re>0, 

is a solution to the following numerical scheme with presumably small forcing terms and zero initial data: 

Eo-=0 (^j,n+a + E(t= 0 E?=—r ^j+i,n+(7 — ^t^j,n+k , j ^ f i R ^ 0 , 


(3.1) 


ej,n — 


Vj,r. 


e-i.o = ■ ■ ■ = Cj^k-i = 0 


ei,o 


0<j<r — 1, n > k , 

j >0. 


The aim of the following two subsections is to quantify the smallness of the source terms in (3.1) in order 
to apply the stability estimate of [GT81]. The smallness of the source terms will yield, up to losing some 
powers of At, a semigroup estimate for (ej,n)j>o,n>o which will eventually give the semigroup estimate 
for the numerical solution {u'j)j>o^n>o- 
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3.1. The case of an incoming velocity. We assume here a > 0 so that no boundary layer arises in 
the solntion to (1.4), (1.6) (Cnum = {0}). The approximate solution merely reads: 

^ / uo{y - a^-) dy, i>0, n>0, 

Jxj 

where we recall that uq has been extended by zero to M“. From the flatness conditions wo(0) = 'Uo(O) = 0, 
we have uq G Ff^(]R). The errors in (3.1) satisfy the following bounds. 

Proposition 3.1. Let us assume a > 0. Under the assumptions of Theorem 1.5 and in the CFL 
regime (1.5), there exists a eonstant C > 0 that is independent of uq and At G (0,1] such that 

(3.2) sup < CAt^ ||no||^2(R+), 

n>k 

r—1 

(3.3) At < C ||uq||^2(r+) . 

n'>k j=0 


Proof. Let us first consider the boundary error terms {r]j^n)- Since uq vanishes on M“, there holds pj^n = 0 
if n > r/(a A). The sum in (3.3) therefore rednces to finitely many terms (and the number of such terms is 
independent of At). We consider some space index j G {0,..., r — 1} and some time index k < n < rjia A), 
and write 


r]j,n = ^ / no(x - ot”)da: = — / / UQ{y)dydx. 

o X j ^ X j 0 


We then apply the Canchy-Schwarz ineqnality and get 


h 




< c 


Tj + l 


px—a 

Uo{yfdydx<CAt\\uo\\l2(^^+) 


Summing the finitely many nonzero error terms, we get (3.3). 

We now deal with the consistency error in the interior domain. Using the consistency assumptions 1.1 
and 1.2, we have: 


^j,n+k . 


2,^0 / tto(x — at"'”^'^) — ito(x — ot”) dx 

f7=0 

\ ^-1 P / r^j+l 

+ ^^/3<7 / uo(x-at”+‘^)dx- / uo{x - 

a=0 l=-r 

1 ^ /-o 

= -— / / Uo{x + y - at'^)dydx 

_(v V X q t/ — Qj (J i 


a=0 

k-1 p 


\ ^ p rxj+i p£Ax 

+ Pa ai / UQ^x + y-a t’^+'^) dy dx 

„-n 0- ^ Jxi Jo 


2^ 

Ax 


<7=0 i=—r 
k 


E - rxj + i rlAx 

aa^a / Uq{x — aa Xy — at"') dy dx 

(7=0 ^ 

\ ^ + i pAx 

+ 2;- Pa ^at / Mo(x + £y-at''+'")d?/dx. 

,7=0 e=-r "'0 



Using the consistency assumptions 1.1 and 1.2 again, we can add the zero quantity 





l=—r 



Uq(x — at") dydx . 
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and get 


^j,n+k — A 


E " r^j+i ^ ^ 

a a^a / / Uq{x + x' — a t^) dx dy dx 

^_n Jxn Jo J—acjXy 


Ax 


(3.4) 


+ 


Ax 


cr=0 
k—1 p 


Xj^i pAx riy—aa At 


EaE<»'/ / / 

_n « _ Jx^ Jo Jo 


Uq {x + x' — at"') dx' dy dx . 


(T=0 i=—r '' 

We now apply successive Cauchy-Schwarz inequalities. In the CFL regime (1.5), we get for instance 

2 


rxj + i pl\x pO 

/ / / Mq(x + x'— at”) dx'dy dx 

J Xj Jo J —a cr \y 


< CAt^ 

< CAP 

< CAP 


rxjj^i pAx pO 


f Xj Jo J —a (7 Xy 

rxj+i pO 


u'q(x + x' — a t")^ dx' dy dx 


pxj+i pu 

/ / Poi^ + x' — a t")^ dx' dx 

Jxj J—acrAt 

/ Uq{x — at")^ dx . 

J x-j—a k At 


The other error term in (3.4) is estimated similarly, and in the end, we can show that there exists a fixed 
integer jo > 0 (that only depends on the CFL number A, a and k) such that 

Vj>r, VuGN, Isj^n+kl"^ a C At / u'o'(x — at”)^ dx. 

The estimate (3.2) follows immediately. □ 


3.2. The case of an outgoing velocity. From now on, we consider the case of an outgoing velocity 
a < 0 for which non-trivial boundary layers appear in the solution to the numerical scheme (1.4), (1.6) 
(Cnum = 1^)- The following Proposition provides error bounds for the source terms in the numerical 
scheme (3.1). 

Proposition 3.2. Under the assumptions of Theorem 1.5 and in the CFL regime (1.5), there exists a 
constant C > 0 that is independent of uq and At G (0,1] such that 


(3.5) 

sup ^ Ax|ej>p < CAP ||u'o'||2 2 (r+) , 

n>2k 

(3.6) 

sup ^Ax|ej>p < CAt ||u'o||^i(]R+), 

k<n<2k—l 
- - 3>r 


r—1 

(3.7) 

< C At Uo . 


n>k j=0 


Proof. We first prove (3.7) and then deal with (3.5) and (3.6). 

Errors at the boundary. We start with the proof of the estimate (3.7). From the definition (2.14), we 
obtain (recall n > k and j = 0,..., r — 1): 


1 r^j+i 1 /■*”'+' 

'nj,n = uf^^ = — / 1 x 0 ( 3 : + |a| t”)dx - — / uo(|a|i)df. 

JXj Jt^ 


With the notation (1.5), the error rjj^n can be written as 

fAt/' 


T r^t/T I 

hj,n = -j^ / ito(xj + |a| t” + |a| s) — ito(|a| t”) ds — — / uo(|fl| + |a| s) — iio(|a| t”) ds 
At Jo At Jq 

fAt/r pxj+\a\s i pi 

Jq UQ{\a\t" + y)dyds- — 


/o ^0 


-I pAt p\a\ s 

' ' UQ{\a\t" + y)dyds. 

'0 Jo 
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Each term in rjj^n is estimated by applying the Cauchy-Schwarz inequality. For instance, we have 

2 


Aa\ s 


Jq Jq 

and similarly, we have 

^ j-M/r j^Xj + \a\s 

If) If) 


Uo{\a\r + y) dyds 


p\a\ At 

<CAt Uo{\a\t^+ yfdy, 

Jo 


no(|a|P + y) dyds 


pr Ax 

<CAt Uo(|a| t"-+ dy . 

Jo 


Summing over the n’s and the finitely many j’s, we derive the bound (3.7). 

Errors in the interior. We decompose the consistency error Sj^n in (3-1) as 

p. = pint I bl 

with self-explanatory notation. The estimate of the interior consistency error follows from the exact 
same arguments as we used in the case of an incoming velocity. The only difference is that, because of 
the sign of a, we do not need to extend uq by zero to and no assumption on the behavior of uq at 0 
is needed to derive the estimate 


(3.8) 


sup ^ < CAt‘^ \\uo\\l 2 ^^+j. 


n>k 
_ 


We now focus on the new consistency error that comes from the boundary layer terms in 


k-l 


-bl 




At 




bl,0 


1 ^ 


fc-i 


bl,l 


\cr=0 

(T = 0 

£=—r 

/ f7=0 

k 

sr ^ bi,o , P 

Y 

k 

k-l 

P 

Yf 

bid , bl,l 

^i,n+fT + Y Y ^j+£,n+a 


+ Y 

cr=0 i=—r 


bl,l 

^i+£,n+(T 


(j=0 (T=0 (T=0 i=—r 

In the case n > k, we use the definition of the boundary layer profile and corrector to simplify 

the latter expression and get^ 

k k—1 p 


-bl 


1 


~-j,n+k ^ 




bl,l _ bl,l 
'j+^,n+(J ^j+£,n 


a=0 


CT=0 £=—r 


The first boundary layer corrector is given in subsection 2.5. In particular, the error can be 

decomposed as a linear combination of sequences of the form 

Wj+l ^ 


At 


'y ^ (Un+(T-|-(T' ^n+cr') ) 


( t '=0 


with a = 0,..., k and i = —r, ■ ■ ■ ,p. Since the sequence (wj) is exponentially decreasing, we have 


n>k 
_ 


n>k ^ ^ 

— cr=0 


sup Y ki.n+fcP ^ sup ^ — 

0^ 

k 

Y^-'( 


C 

< sup > 

At n>k n 

— (T=0 


Y1 l^n+a+a' 

a'=0 


' l^n+a+a' 

tr 

- ^n+a‘ 


.tr 


.tr 


tr> 


a '=0 


We compute 


fAl per At pa'At 


“n+cr+cr' “n+cr' “n+tr ' “n 


'0 JO 


; ea £\t 

/ '^Jo(|o| C -|- |a| Si + |a| S2 + |a| S3) dss ds2 dsi 

Jo 


^Here we use the consistency assumption 1.2. 
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and the Cauchy-Schwarz inequality yields 


,tr 


f(2fc+l) At 


I'^n+fT+o-' ’^n+cr' ''^n+a '^n\ ^ C At 


ito(|a| t" + |a| ds . 


We have thus derived the estimate 


sup Ax\£^]n+k\^ < CAt^ \\uo\\h(m+) ■ 


n>k 

- J>r 


Together with (3.8), this already proves (3.5). 

We turn to the proof of (3.6). It only remains to estimate the norm of (e^).), • • ■, (e^2fe-i) because 
the interior errors (e™fc),. ■ •, (e™2fc-i) already satisfy (3.8), which is not larger than the right hand side in 
(3.6). Let us explain how we derive the estimate for (e^fc)- The remaining terms are similar. The error 
reads 


_ 1 bl,0 , 1 bl,l _ 'Wj tr Wj 


k-l 


M _ 




(T = 0 


cr=0 


so we have 


bl |2 < 


E Ko 


„,tr _ 
^k+a ~ 


u'oiy) dy ds. 


j^r cr=0 

We now use the assumption no(0) = 0 of Theorem 1.5 and get 

^ nAt Ha| (t’^+'^+s) 

Jo Jo 

The Cauchy-Schwarz inequality then gives 

rial i'=+“'+l 

\u^k+a\‘^ < C At Uo(?/)^d?/ < C At^ II“oIIl°°(R+) — C lkc)llL7i(R+) • 

We thus get (3.6) for the sequence {£j^k) and the remaining terms (e^),.,.]^),..., (e^2A:-i) dealt with in 
the same (rather crude) way. □ 

Propositions 3.1 and 3.2 imply the following result which uses GKS type norms. 

Proposition 3.3. Under the assumptions of Theorem 1.5 and in the CFL regime (1.5), there exists a 
constant C > 0 that is independent of uq and At G (0,1], such that for all 'y > 0 there holds 


(3.9) 


(3.10) 


EE At Ax e 


— 2 n 7 At 


-j,n| ^ 


n>k j>r 


<C [1 + -] AU 

7. 


n'>k j=0 


'j,nJ < C At^ 


Proof of Proposition 3.3. The proof of (3.10) is immediate and follows from either (3.3) or (3.7) by using 
7 > 0 (so that the exponential factors in (3.10) are not larger than 1). 

The proof of (3.9) follows from either (3.2) or (3.5)-(3.6). In the incoming case (a > 0), we use (3.2) 
and get 


yy At Axe 


—2n'y At 




< C At^ ||rto| 


n>k j>r 


H^{R+) 

n>k 

< 


'Y 

C 

yi'y At _ y 


Ac 11 Wo 


— ^ II'“o|Ih 2 (r+) 
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which is even better than (3.9). In the outgoing case, we use (3.5)-(3.6) and get 


yy At Axe 


—2 717 Ai 




< C ||uo| 


/J 2 (R+) 


+ C At^ ||uo| 


n>k j>r 


H 2 (R+) 

n>2 k 


p-2n7At 

i>2k 

<c + -^ At^ ||uo||^2(r+) 


□ 


Remark 3.4. If we had not included the boundary layer eorrector in the approximate solution, the right 
hand side in the error estimate (3.9) would have been of the form At ||iJo||^ 2 (]k+) instead of At^ II^o||^ 2 (]r+); 
which would have not been sufficient to derive (1.16) beeause there is a loss of a faetor At in the derivation 
of the estimate (3.12) below. 

3.3. The semigroup estimate. We now prove Theorem 1.5. We apply the main result® of [GT81] which 
states that the numerical scheme (3.1) is strongly stable in the sense of [GKS72]. In other words, there 
exists a constant C > 0 , that is independent of the parameter 7 > 0 , such that there holds: 


(3.11) 


7 




1 + 7 At 


y^AtAxe-^'^^^^le^p + y^ yy Ate 


n>0 j>0 

1 + 7 At 


n>0 ^=0 


—2n'y At \^n\2 

\^j I 


r—1 


< c 


7 


y^ y^AtAxe-2"^^‘|e”|2 + y^ y^Ate-^’^^^^^T^ 


n>k j'>r 


n>k j=0 
+2 11 , I|2 


< C At ||uo||jy 2 (]K+) 


7 At + 1 

7 


1 

1 + - 
7 


1 


where we have used Proposition 3.3 to derive the second inequality in (3.11). We choose 7 = At^, with 
/i G [0,1/3]. We thus derive from (3.11) the bound 

yy Ate"^”^*^'^'' yy Axje^p < CAt^-^^ ||uo|Ih2(r+) • 

n>0 j>0 

In particular, a very crude lower bound for the left hand side gives 


(3.12) 


sup e 

n>0 


-2nAt^+i^ 


y^Axje^P < CAt^-®^||uo| 

j>0 

The semigroup estimate (3.12) yields the bound 

Vn G N, yy Ax < 2 yy Ax + C At^“®^ ||uo||^ 2 (]r+) , 

i>o i>o 

with a constant C that is uniform with respect to all the parameters. We now derive a semigroup estimate 
for the approximate solution In the case of an incoming transport equation (a > 0), we have 

^ / Mx-a e) dx , 

for all j, n G N (recall that uq vanishes on M“). In particular, the Cauchy-Schwarz inequality yields 


y^Axiu^ppp < 


Ujln I ^ IPO 


J>0 


and we get 


VnGN, '^Ax\u]\‘^<2 

j>0 


L 2 (]R+) + C” 


®As a matter of fact, the main result of [GT81] requires more restrictive conditions than Assumption 1.3, but the extension 
of the result of [GT81] to numerical schemes that satisfy Assumption 1.3 was performed in [Goul3]. 
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which gives (1.16). We now consider the case of an outgoing transport equation (a < 0) and derive a 
semigroup estimate for the approximate solution We still have 

< ||uo||i2(K+) , 

j>0 

and we thus focus on the semigroup estimate for the boundary layer profile and corrector. Let us first 
consider the boundary layer profile for which we have 


sup Ax = sup Y P = sup Ax Y 


Wn 


n>0 ^.>0 


n>k 


n>k 

c 

< sup — 

n>k 


j>0 


rto(|a| t) dt 


/i" 


< C ||mo 


We now deal with the first boundary layer corrector for which we have 


sup 

n>0 
- J>0 


Y, |Axu^' 


bl,l|2 

' = sup 


n>k 

- j>o 


Y, — sup C Ax 


n>k 


E 

a '=0 

k 


Ola '^n+a 


E' 

i>o 


< 


CAt sup Y \Y+Y < C ||uo||i2(K+) • 


n>k n 

— (T=0 


As in the incoming case, we have thus derived the bound 

^ Ax< C ||uo|li2(K+) , 

i>o 

and we get (1.16) accordingly. 


4. Example and counterexample 


4.1. A 4 time-step 5 point centered scheme. As a first numerical illustration of the above results in 
the case of an outgoing velocity a = — 1, we consider the following numerical scheme. The time-stepping is 
solved using the 3rd order explicit Adams-Bashforth method, so that assumption 1.2 is satisfied. The space 
discretization of the advection term a dxU is based on a centered five-point approximation supplemented 
with a fourth order stabilizing dissipative term: 


u 


n+l _ n 


= n?-A 


(4.1) 


m _ m—1 ^ m—2 


f? ■■= a 


- 8u]_, + n-_2 + 4n-+, - 6u] + 


12 24 

As we will show with numerical experiments, this scheme displays numerical boundary layers when com¬ 
bined with Dirichlet boundary conditions. Let us compute A: 

A{z) = + 8z- 8z~^ + zY - + 4z-6 + 4z~^ - zY, 

from which we get A(l) = 0 and A'{1) = a and the space discretization satisfies assumption 1.1. Moreover, 
for any 0 G M, one gets 

aY) = —a^(sin(20) — 8sin(0)) — t^(— cos(20) -|- 4cos( 0) — 3), and 51?(A(e*®)) = ^ sin^ ( ^ ) , 


6 


12 


Therefore the only root of A(e*^) in [—7r,7r] is 0 = 0. This ensures that assumption 2.1 is satisfied. 
Figure 4.1 below pictures the closed curve {—AA(e*^), r/ G M| for the choice A = 0.4 (blue curve), together 
with the stability domain of the time integrator (red dashed curve) ; see [HW96, HNW93] for details. Let 
us observe that the stability assumption for the Cauchy problem 1.3 is satisfied. 
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Figure 4.1. Verification of the stability assumption 1.3 for the 3rd order scheme (4.1). 

The numerical test case concerns the following initial condition: 

no(x) = xG[0 ,l], 

together with homogeneous Dirichlet conditions at both left and right boundaries (with no significant 
effect arising from the right boundary due to the incoming situation with a = — 1 and the absence of 
boundary layer at x = 1). We compute the solution on N = 216 uniformly spaced grid cells, until time 
T = 0.5. At this time, the initial bump crosses the left boundary with the highest strength. As expected, 
the numerical solution (u”) develops some boundary layer in the neighborhood of x = 0 due to the 
incompatibility of the homogeneous Dirichlet condition n” = 0, 0<j <r — 1, with the effective trace of 
the solution u‘“*(0'’“, T) = 1. We then observe on Figure 4.2 an oscillating pattern that does not disappear 
as Ax tends to 0. The two roots of A in D \ {0} are real and distinct; one of them equals approximately 
0.0809 and therefore belongs to (0,1), while the second one equals approximately —0.6595 and therefore 
belongs to (—1,0), which gives rise to the oscillations in the boundary layer. 



Figure 4.2. Numerical solution and exact solution at time T=0.5 (216 grid points). 

The main term of the boundary layer expansion is a linear combination of two geometric sequences 
generated by the roots of the equation A(z) = 0 in D \ {0} (see Lemma 2.3). In the present case, we 
obtain numerically zi ~ —0.6595 and Z 2 — 0.0809. The precise boundary layer expansion + 

AxM^*’^(j, T) is depicted with crosses on the left picture of Figure 4.3 for the first 20 grid cells. Notice 
that it depends only on the trace of the solution at the considered time and on the discrete in time 
derivative of this trace, through It fits quite well the difference between the numerical solution and 
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the exact one ix” — T). On the right picture of Figure 4.3 is represented the error in this boundary 

layer expansion T) + AxT)] in the first 50 grid cells. 




Figure 4.3. Boundary layer expansion at time T=0.5 (216 grid points) 


The scheme (4.1) is third order in time and space accurate. We now consider the effective accuracy of 
this scheme for the IB VP problem by computing the £^([0,1]) error at a final given time for successive 
values of 2^ grid points, 5 < M < 12. More precisely, given a time T > 0, we compute the following two 
quantities, where n = Nt is the first integer such that Ay At > T: 


/2" \ 

1/2 


( 2 ^ \ 

J^Ax |u^^-u”*(xj,r)f 

1 . 

and 

^Ax |ix”-M’^PP(xj,r)|^ 


At a first time T = 0.125 at which no significant boundary layer has appeared at x = 0, the convergence 
of both quantities occur with order 3, see Figure 4.4 on the left. For very thin grids, one observes however 
a slight loss of accuracy when computing the usual numerical error. It corresponds to the presence of a 
very small boundary layer that deteriorates the effective order of accuracy. 

At a later time T = 0.4 at which the boundary layer is sufficiently high to affect the convergence, the 
usual numerical error is strongly increased and the apparent order of accuracy is severely damaged: in 
Figure 4.4 on the right, we observe a numerical accuracy of order 0.5 for the usual numerical error, and 
of 1.5 for the error in the boundary layer expansion. 


4.2. The leap-frog scheme. We now consider the usual three-time step leap-frog scheme, with a three 
point stencil in space: 


(4.2) 


u 


n+l 


— U 


n—1 


U 


2At 


a ■ 


j+i 


— u 


1-1 


2 Ax 


= 0 . 


The scheme corresponds to the so-called Nystrom method of order 2 (also called the mid-point formula) 
combined with the center differentiation formula for the space discretization. The corresponding function 
A equals z — and therefore vanishes at —1. Assumption 2.1 is no longer satisfied and Figure 4.5 below 
illustrates that the failure of Assumption 2.1 gives rise to a completely different behavior. Namely, we 
compute the numerical solution for (4.2) with a = — 1 and homogeneous Dirichlet boundary conditions 
at different time levels, for the same kind of bump initial data. As the bump crosses the left boundary, a 
highly oscillatory wave packet emerges from the boundary and propagates with velocity -|-1 towards the 
right. The envelope of this wave packet is exactly the one of the initial condition, see Figure 4.5. The 
latter phenomenon has long been identified of course, see, e. g., [Tre82]. 
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Figure 4.4. Convergence in log/log scale. Solution at time T = 0.125 with no significant 
boundary layer (left) / at time T = 0.4 with an important boundary layer (right). 



Figure 4.5. Leap-frog scheme, solution at time T = 0, T = 0.2, T = 0.5 and T = 1 
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